%% the script red monthly compressed file from the IMAGE array


for nyears = 2011:2013,
    
    for nmonths = 1:12,
        
        date_string = sprintf('%d%02d',nyears,nmonths)
        
        S = urlread(['http://space.fmi.fi/cgi-bin/imagecgi/image-month.cgi?start=' date_string '&yourname=Manoj&institute=UniversityOfColorado&email=manoj.nair@colorado.edu&message=Ionosphere']);
        urlwrite(['http://space.fmi.fi/image/plots/image' date_string '.tar'],['/Users/manojnair/projects/obs_mag_data/image/image' date_string '.tar']);
    end;
end;

%% script to untar and un-gzip the files

for nyears = 2010:2013,
    
    for nmonths = 2:12,
        date_string = sprintf('%d%02d',nyears,nmonths)
        
        untar(['/Users/manojnair/projects/obs_mag_data/image/image' date_string '.tar'],'/Users/manojnair/projects/obs_mag_data/image/dailyfiles/');
    end;
end;

%% now unzip the gzip the IAGA binary files


start_date = datenum(1995,1,1);
end_date = datenum(2012,13,31);


for thisdate = start_date:end_date,
    
    gunzip(['/Users/manojnair/projects/obs_mag_data/image/dailyfiles/image.' datestr(this_date,'yymmdd') '.gz'],'/Users/manojnair/projects/obs_mag_data/image/dailyfiles/');
    
end;

%% now time avrerage the 10 sec files to 60 sec (1 min) files



start_date = datenum(1995,1,1);
end_date = datenum(2012,12,31);


for this_date = start_date:end_date,
    
    system(['/Users/manojnair/projects/obs_mag_data/iaga_average.exe' ...
        ' -t 60 /Users/manojnair/projects/obs_mag_data/image/dailyfiles/image.' datestr(this_date,'yymmdd') ...
        ' > /Users/manojnair/projects/obs_mag_data/image/averageddailyfiles/image_1min_iaga.' datestr(this_date,'yymmdd')]);
    
    fprintf('%s\n',datestr(this_date,'yymmdd'));
    
end;

%% now convert the 1 minute data to MATLAB format

start_date = datenum(1995,1,1);
end_date = datenum(2012,12,31);


for this_date = start_date:end_date,
    
    system(['/Users/manojnair/projects/obs_mag_data/iaga_matlab.exe' ...
        ' /Users/manojnair/projects/obs_mag_data/image/averageddailyfiles/image_1min_iaga.' datestr(this_date,'yymmdd') ...
        ' > /Users/manojnair/projects/obs_mag_data/image/averageddailymatfiles/image_1min_iaga.' datestr(this_date,'yymmdd') '.mat']);
    
    fprintf('%s\n',datestr(this_date,'yymmdd'));
    
end;


%% find out howmany unique stations are there in all the files


start_date = datenum(1995,1,1);
%end_date = datenum(2008,12,31);
end_date = datenum(2012,12,31);

fid = fopen('/Users/manojnair/projects/obs_mag_data/image/station_list.txt','wt');
global_station_list = 'TRO';

for this_date = start_date:end_date,
    
    
    
    % load(['/Users/manojnair/projects/obs_mag_data/image/averageddailymatfiles/image_1min_iaga.' ...
    %     datestr(this_date,'yymmdd') '.mat'],['*' datestr(this_date,'yymmdd')]);
    
    S = whos('-file',['/Users/manojnair/projects/obs_mag_data/image/averageddailymatfiles/image_1min_iaga.' ...
        datestr(this_date,'yymmdd') '.mat'],['*' datestr(this_date,'yymmdd')]);
    
    %     S = whos(['*' datestr(this_date,'yymmdd')]); % Get all the data variables
    %
    array = [S(1:end).size];
    
    
    if (all([S(1:end).bytes]==34608) && all(array(1:2:end)==1442) && all(array(2:2:end)==3))
        
        for i = 1:length(S),
            
            this_stations(i,:) = S(i).name(1:3);
        end;
        
        global_station_list = union(global_station_list, this_stations,'rows');
        
    else,
        fprintf('Error in station %s\n', datestr(this_date,'yymmdd'));
        
    end;
    
    
    %
    %     clear(['*' datestr(this_date,'yymmdd')]);
    %
    %
    fprintf(fid, '%s\n',datestr(this_date,'yymmdd'));
    fprintf(fid, '%s\n',this_stations);
    
    
    clear this_stations;
    
    
    
end;


%% We now have the global station list. Now create a time series for each
% station,

%define global station data from the image array
global_station_list = ...
    [ 'ABK';  'AND';  'BJN';  'DOB';  'DON';  'HAN';  'HOP';  'HOR';  'IVA';...
    'JCK';  'KAR';  'KEV';  'KIL';  'KIR';  'LEK';  'LOZ';  'LYC';  'LYR';...
    'MAS';  'MEK';  'MUO';  'NAL';  'NOR';  'NUR';  'OUJ';  'PEL';  'RVK';...
    'SOD';  'SOL';  'SOR';  'TAR';  'TRO';  'UPS'];
%define a time array
time_array_min = (datenum(1995,1,1,0,0,30): (1/(24*60)): datenum(2012,12,31,23,59,30))';
a= length(time_array_min);



%define start and end dates
start_date = datenum(1995,1,1);
end_date = datenum(2012,12,31);

% fid = fopen('/Users/manojnair/projects/obs_mag_data/image/station_list.txt','wt');
% global_station_list = 'TRO';

for i = 31:length(global_station_list), %for all the stations
    
    x_data = nan([1,a]);
    y_data = nan([1,a]);
    z_data = nan([1,a]);
    
    for this_date = start_date:end_date, % for all the dates
        
        %check whether the station variable availabe from the daily file
        S = whos('-file',['/Users/manojnair/projects/obs_mag_data/image/averageddailymatfiles/image_1min_iaga.' ...
            datestr(this_date,'yymmdd') '.mat'],[global_station_list(i,:) datestr(this_date,'yymmdd')]);
        
        
        if length(S) == 1, %OK data are available. Now load it.
            
            start_index = (this_date - start_date ) * 1440 + 1; % Start index for sotring the array
            
            load(['/Users/manojnair/projects/obs_mag_data/image/averageddailymatfiles/image_1min_iaga.' ...
                datestr(this_date,'yymmdd') '.mat'],[global_station_list(i,:) datestr(this_date,'yymmdd')]);
            
            eval(['data = ' global_station_list(i,:) datestr(this_date,'yymmdd') ';']);
            
            x_data(start_index:start_index + 1440 - 1) = data(3:end,1);
            y_data(start_index:start_index + 1440 - 1) = data(3:end,2);
            z_data(start_index:start_index + 1440 - 1) = data(3:end,3);
            
            clear([global_station_list(i,:) datestr(this_date,'yymmdd')]);
            
        end;
        
    end;
    
    subplot(311);
    
    plot(time_array_min,x_data);
    datetick('x');
    ylabel('X (nT)');
    title(global_station_list(i,:));
    subplot(312);
    plot(time_array_min,y_data);
    datetick('x');
    ylabel('Y (nT)');
    subplot(313);
    plot(time_array_min,z_data);
    datetick('x');
    ylabel('Y (nT)');
    
    saveas(gcf,['/Users/manojnair/projects/obs_mag_data/image/plots/' global_station_list(i,:)],'png');
    eval(['save /Users/manojnair/projects/obs_mag_data/image/station_matlab_files/' global_station_list(i,:) ' x_data y_data z_data']);
    
    close all;
end;


%% Now create netCDF files for these stations

[master_a,master_b] = xlsread('/Users/manojnair/projects/obs_mag_data/image/work/IMAGE_ARRAY.xlsx');
master_b(1,:) = [];

% get the total stations

localdirectory = '/Users/manojnair/projects/obs_mag_data/image/Absolute_netCDF/';

time_array_min = (datenum(1995,1,1,0,0,30): (1/(24*60)): datenum(2012,12,31,23,59,30))';
a = length(time_array_min);

for i = 1:33,
    
    
    metadatacellarray = master_b(i ,:);
    lat = master_a(i ,1);
    long = master_a(i ,2);
    alt = master_a(i ,3);
    fprintf('Creating %s \n',cell2mat(master_b(i,1)));
    [ncid, X_ID, Y_ID, Z_ID] = Create_Geomag_netCDF(localdirectory, lat, long, alt, metadatacellarray);
    eval(['load /Users/manojnair/projects/obs_mag_data/image/work/station_matlab_files/' cell2mat(master_b(i ,1)) '.mat' ' x_data y_data z_data']);
    x_data(isnan(x_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
    y_data(isnan(y_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
    
    z_data(isnan(z_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
    
    %multiplied by 10 to preserve the decimal points
    
    netcdf.putVar(ncid, X_ID, 0, a,  int32 (x_data .* 10) );
    
    netcdf.putVar(ncid, Y_ID, 0, a,  int32 (y_data .* 10) );
    
    netcdf.putVar(ncid, Z_ID, 0, a,  int32 (z_data .* 10) );
    
    netcdf.close(ncid);
end;


%% Now fit spline, calculate IGRF, remove the secular variation and create the residual
% netCDF files


[master_a,master_b] = xlsread('/Users/manojnair/projects/obs_mag_data/image/work/IMAGE_ARRAY.xlsx');
master_b(1,:) = [];

% get the total stations

localdirectory = '/Users/manojnair/projects/obs_mag_data/image/Secular_Removed_netCDF/';

time_array_min = (datenum(1995,1,1,0,0,30): (1/(24*60)): datenum(2012,12,31,23,59,30))';
a = length(time_array_min);

for i = 18:18,
    
    
    metadatacellarray = master_b(i ,:);
    lat = master_a(i ,1);
    long = master_a(i ,2);
    alt = master_a(i ,3);
    fprintf('Creating %s \n',cell2mat(master_b(i,1)));
    [ncid, X_ID, Y_ID, Z_ID] = Create_Geomag_netCDF(localdirectory, lat, long, alt, metadatacellarray);
    eval(['load /Users/manojnair/projects/obs_mag_data/image/work/station_matlab_files/' cell2mat(master_b(i ,1)) '.mat' ' x_data y_data z_data']);
    %calculate igrf
    
    % For some stations selectiely delete the bad data
    
    if cell2mat(master_b(i,1)) == 'LYR',
        L = time_array_min >= datenum(2006,1,1) & time_array_min <= datenum(2007,8,1) | ...
            time_array_min >= datenum(2008,8,1);
        x_data(L) = NaN;
        y_data(L) = NaN;
        z_data(L) = NaN;
    elseif cell2mat(master_b(i,1)) == 'MUO',
        L = time_array_min >= datenum(1995,1,1) & time_array_min <= datenum(2005,1,1);
        x_data(L) = NaN;
        y_data(L) = NaN;
        z_data(L) = NaN;
    elseif cell2mat(master_b(i,1)) == 'PEL',
        L = time_array_min >= datenum(1995,1,1) & time_array_min <= datenum(1997,1,1);
        x_data(L) = NaN;
        y_data(L) = NaN;
        z_data(L) = NaN;
    elseif cell2mat(master_b(i,1)) == 'KIL',
        L = time_array_min >= datenum(2009,1,1) & time_array_min <= datenum(2013,1,1);
        x_data(L) = NaN;
        y_data(L) = NaN;
        z_data(L) = NaN;
    elseif cell2mat(master_b(i,1)) == 'IVA',
        L = time_array_min >= datenum(2012,1,1) & time_array_min <= datenum(2013,1,1);
        x_data(L) = NaN;
        y_data(L) = NaN;
        z_data(L) = NaN;
    end;
    
    
    icount = 1;
    for decyear = 1995:0.5:2013,
        mag_igrf(icount,:) = igrf11magm(0,master_a(i,1),master_a(i,2),decyear);
        icount = icount + 1;
    end;
    
    %fit spline
    
    L = isnan(x_data);
    data_array = x_data(~L);
    time_array = time_array_min(~L);
    b1 = min(time_array):365:max(time_array)+10;
    if length(b1) > 1,
        sp=spline(b1,data_array(1:60:end)/spline(b1,eye(length(b1)),time_array(1:60:end)));
        v=ppval(time_array_min,sp);
    else,%data length <=1 year
        sp=robustfit(time_array,data_array);
        v= sp(1) + sp(2) * time_array_min;
    end;
    subplot(311);
    plot(time_array_min,x_data);
    hold on;
    plot(time_array_min(~L),v(~L),'r');
    plot(datenum(1995:0.5:2013,1,1),mag_igrf(:,1),'c.');
    ylabel('X');
    axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
    datetick('x','keeplimits');
    title(cell2mat(master_b(i,1)));
    % write sec removed z values
    x_data = x_data - v';
    
    L = isnan(y_data);
    data_array = y_data(~L);
    time_array = time_array_min(~L);
    b1 = min(time_array):365:max(time_array)+10;
    if length(b1) > 1,
        sp=spline(b1,data_array(1:60:end)/spline(b1,eye(length(b1)),time_array(1:60:end)));
        v=ppval(time_array_min,sp);
    else,%data length <=1 year
        sp=robustfit(time_array,data_array);
        v= sp(1) + sp(2) * time_array_min;
    end;
    subplot(312);
    plot(time_array_min,y_data);
    hold on;
    plot(time_array_min(~L),v(~L),'r');
    plot(datenum(1995:0.5:2013,1,1),mag_igrf(:,2),'c.');
    ylabel('Y');
    axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
    datetick('x','keeplimits');
    
    % write sec removed z values
    y_data = y_data - v';
    
    
    L = isnan(z_data);
    data_array = z_data(~L);
    time_array = time_array_min(~L);
    b1 = min(time_array):365:max(time_array)+10;
    if length(b1) > 1,
        sp=spline(b1,data_array(1:60:end)/spline(b1,eye(length(b1)),time_array(1:60:end)));
        v=ppval(time_array_min,sp);
    else,%data length <=1 year
        sp=robustfit(time_array,data_array);
        v= sp(1) + sp(2) * time_array_min;
    end;
    subplot(313);
    plot(time_array_min,z_data);
    hold on;
    plot(time_array_min(~L),v(~L),'r');
    plot(datenum(1995:0.5:2013,1,1),mag_igrf(:,3),'c.');
    ylabel('Z');
    axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
    datetick('x','keeplimits');
    
    % write sec removed z values
    z_data = z_data - v';
    
    
    
    
    x_data(isnan(x_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
    y_data(isnan(y_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
    
    z_data(isnan(z_data)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
    
    %multiplied by 10 to preserve the decimal points
    
    netcdf.putVar(ncid, X_ID, 0, a,  int32 (x_data .* 10) );
    netcdf.putVar(ncid, Y_ID, 0, a,  int32 (y_data .* 10) );
    netcdf.putVar(ncid, Z_ID, 0, a,  int32 (z_data .* 10) );
    
    netcdf.close(ncid);
    
    saveas(gcf,['/Users/manojnair/projects/obs_mag_data/image/plots/' cell2mat(master_b(i,1)) '_Spline_IGRF_fit'],'png');
    close all;
end;


%% Now test read some netCDF files

time_array_min = (datenum(1995,1,1,0,0,30): (1/(24*60)): datenum(2012,12,31,23,59,30))';
a = length(time_array_min);

%ncfname = '/Users/manojnair/projects/obs_mag_data/image/Absolute_netCDF/LYR_1995_2012_Absolute.nc';
ncfname = '/Users/manojnair/projects/obs_mag_data/image/Secular_Removed_netCDF/TRO_1995_2012_Secular_Removed.nc';


[x_data, y_data, z_data , dummy, dummy, dummy, dummy] = read_geomag_netcdf(ncfname, 0, a, 0);

subplot(311);
plot(time_array_min, x_data);
axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
datetick('x','keeplimits');
subplot(312);
plot(time_array_min, y_data);
axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
datetick('x','keeplimits');
subplot(313);
plot(time_array_min, z_data);
axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
datetick('x','keeplimits');

        

